Lattice Boltzmann Algorithm for three-dimensional liquid crystal 

hydrodynamics 

C. Denniston 1 , D. Marenduzzo 2 , E. Orlandini 3 , and J.M. Yeomans 2 
1 Department of Applied Mathematics, 
The University of Western Ontario, London, Ontario N6A 5B8, Canada 
^Department of Physics, Theoretical Physics, 
1 Keble Road, Oxford, 0X1 3NP, England 
3 INFM, Dipartimento di Fisica, Universita' di Padova, Via Marzolo 8, 35131 Padova, Italy 

Abstract 

We describe a lattice Boltzmann algorithm to simulate liquid crystal hydrodynamics in three dimensions. 
The equations of motion are written in terms of a tensor order parameter. This allows both the isotropic and 
the nematic phases to be considered. Backflow effects and the hydrodynamics of topological defects are 
naturally included in the simulations, as are viscoelastic effects such as shear-thinning and shear-banding. 
We describe the implementation of velocity boundary conditions and show that the algorithm can be used to 
describe optical bounce in twisted nematic devices and secondary flow in sheared nematics with an imposed 
twist. 
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I. INTRODUCTION 



Liquid crystals are fluids, typically comprising long thin molecules. Subtle energy-entropy bal- 
ances can cause the molecules to align to form a variety of ordered states. For example, in nematic 
liquid crystals, the molecules tend to align parallel giving a state with long-range orientational or- 
der. Liquid crystals exhibit both an elastic and a viscous response to an external stress. Coupling 
between the director and velocity fields leads to strongly non-Newtonian flow behaviour such as 
shear-banding and molecular tumbling under an applied shear. 

Given both the rich rheological behaviour of liquid crystals and their importance in optical de- 
vices there is need to develop a robust numerical method to allow us to explore their dynamics. 
The hydrodynamic equations of motion are complex and, although much interesting analytic work 
has been carried out (see de Gennes & Prost 1993, Beris & Edwards 1994, and the recent review 
Rey & Denn 2002), there are many interesting questions that cannot be answered analytically. 
Therefore, in this paper, we describe a lattice Boltzmann approach to liquid crystal hydrodynam- 
ics in three dimensions. This generalizes a two dimensional algorithm which was proposed by 
Denniston et al. (2000). For a review of the lattice Boltzmann algorithm see Chen & Doolen 
(1998), S. Succi (2000), Wolf-Gladrow (2000), R. Benzi et al. (1992). 

We first summarize the Beris-Edwards equations of motion for nematic liquid crystals (Beris 
et al., 1990,1994). The lattice Boltzmann algorithm is described in some detail. We pay specific 
attention to boundary conditions which can be more cumbersome in three dimensions than in two. 
We then present two examples, switching in a twisted nematic device and shear flow in a twisted 
nematic, to illustrate the applicability of the method. 

H. MODELLING LIQUID CRYSTAL HYDRODYNAMICS 
A. Equations of motion 

Liquid crystals are described in terms of a local tensor order parameter Q, that is related to 
the direction of individual molecules m by Q a/3 = (m a rhfs — \8 a p), where the angular brackets 
denote a coarse-grained average. Greek indices will be used to represent Cartesian components 
of vectors and tensors, and the usual summation over repeated indices will be assumed. Q is a 
traceless symmetric tensor. Its largest eigenvalue, |g, < q < 1, describes the magnitude of order 
along its principle eigenvector n, referred to as the director. 
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The order parameter evolves according to the equation (Beris et al, 1990,1994) 

(d t + u- V)Q-S(W,Q) = TH (HI) 

where T is a collective rotational diffusion constant and 

S(W,Q) = (£D + n)(Q + I/3) + (Q + I/3)(£D-fi) 

-2f(Q + I/3)Tr(QW) (11.2) 

where D = (W + W T )/2 and fl = (W — W T )/2 are the symmetric part and the anti-symmetric 
part respectively of the velocity gradient tensor W a p = dpu a . The mixture of upper and lower 
convective derivatives is governed by the constant £, which depends on the molecular details of a 
given liquid crystal. 

The term on the right-hand side of Eq. dII.ll> describes the relaxation of the order parameter 
towards the minimum of the free energy. The driving motion is provided by the molecular field 
H, which is related to the variational derivative of the free energy T by 

H=-g + (I/S)7vg 0L3) 

The symmetry and zero trace of Q (and H) is exploited for simplification. 
The fluid has density p and obeys both the continuity equation 

d tP + d aP u a = o (n.4) 

and the Navier-Stokes equation 

pTf 

pd t u a + pupdpUa = d/3T al3 + dpa a p + -^-{dp{{8 a p - 3<9 p P )<9 7 u 7 + d a up + dpu a ). (11.5) 

The form of this equation is similar to that for a simple fluid. However, the details of the stress 
tensor reflect the additional complications of liquid crystal hydrodynamics. There is a symmetric 
contribution 
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&ap — —Po$al3 — ^H Q1 (Qryp + -S^p) — £,(Q a ~f + -5 ay )H 7 p 



+2£(Q Q/ 3 + -5 a p)Q le H^ e - dpQ lv — - (11.6) 
and an antisymmetric contribution 

T ap — Qa-yH^p — H 0ll Q 1 p. (II.7) 



The hydrostatic pressure P is defined below, while rj = ^f- is an isotropic viscosity (77 is related 
to the details of the lattice Boltzmann algorithm, see Eq. III. 1 6b . If more elastic constants are 
introduced in the theory, r af 3 also contains the antisymmetric contribution from dpQ yu sq S q . 

The lattice Boltzmann algorithm described in the next subsection can be used for any model of 
the above form. For the examples in the last section we will use a specific model determined by a 
Landau-de Gennes free energy (de Gennes & Prost, 1993; Doi & Edwards, 1989) 

T = \d\ {^(1 - l)Ql, - ^fQ a ,Q, y Q ja + ^(Ql,) 2 + ^d a Q, x y) (II.8) 

where A is a constant and k denotes the elastic constant of the liquid crystal. We shall work within 
the one elastic constant approximation. Although it is not hard to include more general elastic 
terms this simplification will not affect the qualitative behaviour. The free energy (III. 81) describes 
a first order transition from the isotropic to the nematic phase, controlled by the parameter 7. The 
hydrostatic pressure P is taken to be 

P = pT - |(VQ) 2 , (119) 
and is constant in the simulations to a very good approximation. 



6. Lattice Boltzmann algorithm 

Usually lattice Boltzmann algorithms, describing the Navier-Stokes equations of a simple fluid, 
are defined in terms of a single set of partial distribution functions, the scalars fi(x), that sum on 
each lattice site x to give the density (Chen & Doolen 1998). For liquid crystal hydrodynamics, 
this must be supplemented by a second set, the symmetric traceless tensors Gi(x), that are related 
to the tensor order parameter Q. Each Gj is associated with a lattice vector e*j. We choose a 
15-velocity model on the cubic lattice with lattice vectors: 

ef = (0,0,0) (11.10) 
ef = (±1,0,0),(0,±1,0),(0,0,±1) (11.11) 
ef = (±1,±1,±1). (11.12) 

The indices, i, are ordered so that i = corresponds to ef \ i — 1, • • • 6 correspond to the e*p set 
and i = 7, 14 to the ef 1 set, as illustrated in Fig.[l] 
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FIG. 1: Lattice geometry and lattice vectors for the three dimensional lattice Boltzmann model. 



Physical variables are defined as moments of the distribution function: 

P = Y,fi> P u <* = fi e i°» Q = J2 G i (11.13) 

i i i 

The distribution functions evolve in a time step At according to 

f i (x + e l At,t + At) - fi(x,t) = -1 [C fi (x,t, {f}) + C fi (x + e i At,t + At, {/*})] (11.14) 

Gt(x + e { At, t + At) - Gi(x, t) = -± [C Gi (x, t, {Gj}) + C Gi (x + %At, t + At, {G*})] (11.15) 

This represents free streaming with velocity e*j followed by a collision step which allows the distri- 
bution to relax towards equilibrium. /* and G* are first order approximations to fi(x + e^At, t + 
At) and Gi(x + e^At, t + At) respectively. They are obtained by using AtCti(x, t, {fi}) on the 
right hand side of Eq. dll.141) and a similar substitution in Eq. (III.15I) . Discretizing in this way, 
which is similar to a predictor-corrector scheme, has the advantages that lattice viscosity terms are 
eliminated to second order and that the stability of the scheme is improved. 

The collision operators are taken to have the form of a single relaxation time Boltzmann equa- 
tion, together with a forcing term 

C fi (x,t,{fi}) = --(f i (x,t)- tf q (x,t, {fi})) + Pi (x,t, {fi}), (11.16) 
T f 

C Gl (x,t,{Gi}) = -—(Gi(x,t) - G e i q (x,t,{Gi})) + Mi(x,t,{Gi}). (11.17) 

TG 

The form of the equations of motion and thermodynamic equilibrium follow from the choice 
of the moments of the equilibrium distributions fl q and G eq and the driving terms p { and Mj. f? 9 
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is constrained by 

E fi 9 = P' Y fi 9 ^ = P U a, Y fi 9e ia e if3 = ~^aP + PU a Uf3 (11.18) 

i i i 

where the zeroth and first moments are chosen to impose conservation of mass and momentum. 
The second moment of f eq controls the symmetric part of the stress tensor, whereas the moments 
of Pi 

YPi = ®' YVi e ™ = d/3 T <x/3, YP* e i<* e *P = Q (IL19) 

i i i 

impose the antisymmetric part of the stress tensor. For the equilibrium of the order parameter 
distribution we choose 

Y G i 9 = Q) Y G i 9e i<* = Q M "> Y G 7 e ia e iP = QUaUp. (11.20) 

i i i 

This ensures that the order parameter is convected with the flow. Finally the evolution of the order 
parameter is most conveniently modeled by choosing 

Y Mi = TH(Q) + S(W, Q) = H, Y M * e - = (E M *) u * ( IL21 ) 

i i i 

which ensures that the fluid minimises its free energy at equilibrium. Conditions (III.18l H lH.21l) 
are satisfied, as is usual in lattice Boltzmann schemes, by writing the equilibrium distribution 
functions and forcing terms as polynomial expansions in the velocity 

fi 9 = A s + B s u a e ia + C s u 2 + D s u a upe ia eip + E sa(3 e ia eip, 
Gf = 3 S + K s u a e ia + ~L s u 2 + ~N s u a Uf 3 e ia e i/3} 

Pi = TgdpTapeia, 

Mi = R s + S s u a e ia , (11.22) 

where s G {0, 1, 2} identifies separate coefficients for the velocity vectors ef\ The coefficients 
are determined by evaluating the constraints Eqs. dll. 1 81) — (III.2 1 ft and matching terms. In doing 
this, we have made use of the following symmetry relations for the lattice vectors, Eq. dII.12ft . 

E e fa=E e 2 ) = a = 1,2, 3 (11.23) 

i=l i=7 

Y e£?eg ) = 26*0 Y ^tp = ^ ( n - 24 ) 

i=l t=7 

E e ia e if3 e \r) = E e itt e iB e in = (11.25) 
i=l i=7 



E eSegMJei? = 2 M/*i** ( n - 26 ) 
i=i 

E e$e$e$e® = 8A aPvC - 16<^Ac- (H.27) 



where A a @nc — 5^5^ + 8 art 8p£ + 5 Q c<W 

The expansion coefficients are then determined to be: 



A 2 = ^(TrP/3), A 1 = A 2 , A = p-14A 2 , 

£? 2 = p/24, B l = 8B 2 , 

C 2 = ~7^i Ci = 2C 2 , C = — -p, 

D 2 = ■?-, D 1 = 8D 2 
16 

E 2a /3 = — (— 0"a/3 ~ TrP/35 a p), E la p = 8E 2af3 , 

Jo = Q, K 2 = Q/24, Kx = 8K 2 , 

Q 2Q 

L2 = — Li = 2L 2 , L = — — , 

N 2 = -§, N 1= 8N 2 
lb 

R 2 = H/15, Ri = Ro = R 2 

S2 = Si = 8S 2 , 

T 2 = l/24, T 1 = 8T 2 . (11.28) 



C. Velocity boundary conditions 

To illustrate the implementation of the boundary conditions for a system sheared along y con- 
sider a wall at z — 0. No flux across the wall at z = implies ^ /^e^ = and hence 

h + h + h + h + ho = h + hi + A2 + fn + U (n.29) 

For no slip along the re-direction X)» /ie, x = or 

h + fr + ho + hi + /i4 = /a + /s + / 9 + /12 + hs- (H.30) 
For fixed velocity u* along the y direction J2i h e iy — P u * y or 

h + h + fs + hi + h2 -h-h- ho - fn - /i4 = pu* r (11.31) 
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For the wall at z = there are 5 unknown distributions fn,f7, fs, f9 and fiQ. To determine 
these, two further constraints are required in addition to the relations (III.29I ). (III.30I) . and (III .31b . 
Symmetry suggests 

h~k = fio - U (H.32) 
The conservation of mass is easily obtained by choosing 

fs = U (11.33) 

Other choices can be made for these two further constraints. With the present one, equations 
(III.29I) - (III.33I) can now be solved to give: 

/5 = fe, 

f7 = !(-/l-/2 + /3 + /4-/u + A2 + 3/l3 + /l4 + K)> 

fs = i(/l-/2-/ 3 + /4 + /ll-/l2 + /l3 + 3/i4 + K)' 

h = J(/l + /2-/ 3 -/4 + 3/ii + / 12 -/i3 + /l4-K)» 

/io = \ (-fi + h + h - h + fu + 3/i 2 + /is - fu - K) • ( IL34 ) 

IE. NUMERICAL RESULTS 

Our main purpose here is to present in detail a numerical algorithm for simulating liquid crys- 
tal hydrodynamics in three dimensions. Consequently, we restrict ourselves to presenting a few 
simple examples, aimed at checking the approach. Further numerical applications are listed in the 
concluding section and will be reported elsewhere. 

To apply the algorithm to systems of practical relevance, we need to map the simulation pa- 
rameters listed in Section (2b) onto physical numbers which can be compared to real experiments. 
This task is not trivial and proceeds as follows. First, we have to fix a physical value for the energy 
scale, which in simulations is given by pT. Second, the size of the simulation lattice along z, 
L z , is matched to the liquid crystal cell thickness. Third, the parameter Y is related to the magni- 
tude of order and the physical value of the Miesowicz viscosity 71 (see De Gennes & Prost 1993) 
via the formula Y = 2q 2 /'j 1 (Denniston et al. 2001). This fixes the time step in the simulation. 
Once the energy, length and time scales are specified, all simulation parameters can be related 
straightforwardly to experimental quantities. 
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As specific examples of applications of the three-dimensional algorithm we propose two hydro- 
dynamic effects in two commonly used twisted nematic cells. The simplest example of such cells 
is obtained when at the two surfaces (perpendicular to the z-axis, say) the director is anchored and 
the anchoring is homogeneous, i.e. the directors lie in the surface plane. If the two anchoring di- 
rections at z = and z = L z are mutually perpendicular, the free energy minimum in the absence 
of other factors is a twisted director field along the z direction. This is the geometry we consider 
(with the two pinning directions making an angle of ±7r/4 with the x axis to maximise the system 
symmetry). 

In devices, this cell is switched on by means of an electric field perpendicular to the surfaces. 
Here we shall study the dynamics when the cell is switched off and relaxes back to the twisted 
state. Therefore we initialize the director field to lie along the z direction in the whole cell apart 
from two small regions (8 lattice sites in our simulations) near the top and bottom boundaries to 
allow a smooth relaxation between surface and bulk order. The numerical parameters used for 
this calculation were as follows: k = 0.05, 7 = 3.5, L z — 91, V — 0.33775, £ = 0.8 and 
Tf = 0.56. By means of the procedure outlined above, the simulation system can be mapped onto 
a physical device of thickness 4.5 /im, with all three Frank elastic constants equal to ~ 9.3 pN, 
while T = 0.58 Poise -1 and the Miesowicz coefficient 7! = 1.3 Poise. 

In Fig. 12 we show how the tilt, 9, of the director in the midplane (z = L z /2) changes as a 
function of time immediately after the device is switched off. (The tilt angle is the angle with 
respect to the xy plane). In the figure we compare the results obtained when backflow is neglected 
and when it is included. It can be seen that backflow has a significant effect on the switching, 
namely it causes the non-monotonicity in the relaxation dynamics of the mid-plane tilt angle. The 
flow in the device causes the director to start its relaxation in the 'wrong' direction. As the system 
is symmetric the mid-plane tilt has to pass through 90° again during the relaxation to the twisted 
state. This is a known effect which is called of optical bounce. Note also that with backflow the 
relaxation is faster. 

Optical bounce has only recently been observed directly in an experiment (Smith et al. 2002, 
Ruan et al. 2002). Our curve compares qualitatively well with the one observed by Ruan et al., 
who use a cell of comparable thickness and 71 value, but with a slightly larger value of the elastic 
constants. In the experiments, the maximum in the tilt angle occurred for a value of 9 ~ 96°, the 
backflow effect lasted ~ 2 ms, and after 10 ms the angle had relaxed down to ~ 60°. Since the 
details of the dynamics depend on the width of the region which is initially oriented along z and 
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on the temperature, we are satisfied with semiquantitative agreement here. 

The optical bounce can be seen within Leslie-Ericksen theory (see e.g. Kelly et al. 1999). 
However with the Beris-Edwards equations solved here we can also investigate the dependence of 
the dynamics on temperature. The bounce gets deeper as the effective temperature is decreased 
for constant k, or if k is increased at constant 7. Interestingly, as the temperature increases and 
approaches the isotropic-to-nematic transition point, the waiting time needed for the system to start 
the bounce increases. This happens because the magnitude of the order near the surface, where the 
initial director deformations are confined, drops significantly allowing two defects to transiently 
form. These defects block the system evolution and slow down the start-up of the dynamics. 
Furthermore, we can assess the effect of tuning the strength of anchoring and the magnitude of 
surface order on the bounce magnitude: while the former has a negligible effect on the result, the 
latter can have a bigger impact, i.e. if the order within the surface is slightly different from the 
typical order within the bulk, the bounce magnitude can change significantly. 

As a second example, we consider another twisted geometry, in which the director is anchored 
parallel on the two boundaries, and a n twist is introduced in the z direction. We consider a shear 
experiment in which the bounding plates move at constant velocities ±u y . Because the director 
does not lie in the shear plane, (the yz plane), the hydrodynamics is complex. In particular, 
experiments and approximate modelling (using the Leslie-Ericksen equations) predict that even 
though the system is sheared along y, a secondary flow appears along the x direction (De Gennes 

6 Prost 1993). 

In Figure |3j we show that the lattice Boltzmann simulations can reproduce this three- 
dimensional feature well. In the same figure, we also show how the director field components n x , 
n y and n z vary across the cell. These results were obtained for simulation parameters k = 0.05, 

7 = 3.5, L z = 51, T = 0.33775, £ = 0.85 and r/ = 0.56, which can be mapped this time onto a 
physical device of thickness 2.25 fxm, with all three Frank elastic constants equal to ~ 3.3 pN. 

IV. DISCUSSION 

In this paper we have described in detail a lattice Boltzmann algorithm to simulate liquid crystal 
hydrodynamics in three dimensions. We have shown that it can, for values of the simulation 
parameters which capture physical nematic properties and device dimensions, reproduce optical 
bounce in a twisted nematic device and secondary flow in a twisted nematic under shear. Future 
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FIG. 2: Evolution of the tilt of the director field at the centre of a twisted nematic cell after the field is 
switched off. Notice that 9 and it — 6 actually represent the same state since the director field has no 'head' 
or 'tail'. 
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FIG. 3: Velocity field (left) and director field (right) as a function of the position across the cell for a tt 
twisted nematic liquid crystal under shear. Note the secondary flow along x in the left figure. The velocity 
profile in the left plot is scaled by the shearing velocity, which was ~ 0.8 mm/s. 

applications will include studying backflow in more complicated device geometries such a pi-cell 
(Jung et al. 2003) and 4-domain twisted nematic cells. It will also be of interest to consider the flow 
of liquid crystals in narrow channels as, with the advent of microchannel technology, quantitative 
experiments may soon become possible. Here, surface anchoring and elastic deformations of 
the director field are expected to heavily impact on the rheological properties (Marenduzzo et al. 
2003). Other areas where this approach will be useful is in the study of shear banding under 
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both constant stress and constant strain and in considering the flow of cholesterics, nematic liquid 
crystals with a director field that form a helix. 
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